Variational integrators, the Newmark scheme, and dissipative systems

Variational methods are a class of symplectic-momentum integrators for ODEs. Using these schemes, it is shown that the classical Newmark algorithm is structure preserving in a non-obvious way, thus explaining the observed numerical behavior. Modifications to variational methods to include forcing and dissipation are also proposed, extending the advantages of structure preserving integrators to non-conservative systems. 1 Variational integrators Variational integrators can be derived from two points of view, either as generating function methods or from a discrete Lagrangian mechanics. See Kane, Marsden, Ortiz and West [2000] and Wendlandt and Marsden [1997] for details. Such integrators are automatically symplectic and preserve all momenta. In addition, the symplectic nature of such schemes implies that they will typically also have excellent energy behavior, in the sense of bounded long-time energy error. A particular example of a variational integrator, which is useful in the analysis of the Newmark scheme, is given by the implicit scheme: qk+1 = qk + (∆t)q̇k + (∆t)2 2 ak q̇k+1 = q̇k + (∆t)ak (1) ak = M −1(−∇V ((1 − α)qk + αqk+1)) which updates positions and velocities (qk, q̇k) at time k to values at time k + 1. The mass matrix is given by M and the potential function by V . 2 The Newmark algorithm The Newmark scheme, originally introduced by Newmark [1959], is a classical time-stepping algorithm popular in structural mechanics codes. The method updates points (qk, q̇k) to (qk+1, q̇k+1) by the implicit relations qk+1 = qk + (∆t)q̇k + (∆t)2 2 [(1 − 2β)ak + 2βak+1] q̇k+1 = q̇k + (∆t) [(1 − γ)ak + γak+1] (2) ak = M−1(−∇V (qk)). This algorithm is very similar to the variational scheme specified by (1). The only difference is that the variational scheme evaluates the accelerations at averaged points, while Newmark averages the accelerations. Whereas the variational scheme (1) is symplectic and momentum preserving, it can be readily shown that Newmark, for most values of γ and β, does not preserve the canonical symplectic form and does not preserve momenta. This is clearly pointed out in Simo, Tarnow and Wong [1992], where the authors state: What may seem surprising is that all of the implicit members of the Newmark family, perhaps the most widely used time-stepping algorithms in nonlinear structural dynamics, are not designed to conserve energy and also fail to conserve momentum. Among the explicit members, only the central difference method preserves momentum. Despite this fact, numerical simulations performed with the Newmark code show some remarkable features. From an analysis on linear systems, it is clear that Newmark with γ < 1/2 will dissipate energy, and with γ > 1/2 will produce energy growth. With γ = 1/2 and arbitrary β, however, the energy can be seen to merely oscillate. This is clear from runs such as that in Figure 1(a), which is a simulation of a two degree-of-freedom non-linear oscillator with the Newmark method. The bounded energy error which is observed in such simulations is typical of symplectic integrators, leading one to suspect that Newmark must be structure preserving in some sense. A precise statement of the geometric nature of Newmark is provided by the following theorem, taken from Kane, Marsden, Ortiz and West [2000]. Theorem 1. The Newmark method with γ = 1/2 and any β is the discrete Euler-Lagrange equations for an implicitly defined discrete Lagrangian Lβd . This result establishes that Newmark is actually variational and symplectic, but not in an obvious sense. There are two main ways to interpret this fact. Firstly, one can use this result to prove the existence of momentum functions and a noncanonical symplectic form which will be exactly conserved by the Newmark scheme. Such functions will, however, only be given implicitly and are likely to be extremely unwieldy. Secondly, there exist exactly conserved quantities which are only slight perturbations of the natural momenta and symplectic form. These can be used to prove that the standard quantities will have bounded error for arbitrarily long times. This is similar in spirit to the analysis of energy error often performed for symplectic integrators. 0 10 20 30 40 50 60 70 80 90 100 0.18 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34